General-Relativistic Resistive Magnetohydrodynamics in three dimensions: formulation and tests 
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We present a new numerical implementation of the general-relativistic resistive magnetohydrodynamics 
(MHD) equations within the Whisky code. The numerical method adopted exploits the properties of Implicit- 
Explicit Runge-Kutta numerical schemes to treat the stiff terms that appear in the equations for small electrical 
conductivities. Using tests in one, two, and three dimensions, we show that our implementation is robust and 
recovers the ideal-MHD limit in regimes of very high conductivity. Moreover, the results illustrate that the code 
is capable of describing physical setups in all ranges of conductivities. In addition to tests in flat spacetime, we 
report simulations of magnetized nonrotating relativistic stars, both in the Cowling approximation and in dy- 
namical spacetimes. Finally, because of its astrophysical relevance and because it provides a severe testbed for 
general-relativistic codes with dynamical electromagnetic fields, we study the collapse of a nonrotating star to 
a black hole. We show that also in this case our results are in very good agreement with the perturbative studies 
of the dynamics of electromagnetic fields in a Schwarzschild background and provide an accurate estimate of 
the electromagnetic efficiency of this process. 
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I. INTRODUCTION 

Magnetic fields play an important role in several astrophys- 
ical scenarios, many of which involve also the presence of 
compact objects such as neutron stars (NSs) and black holes 
(BHs), whose accurate description requires the numerical so- 
lution of the equations of general relativistic magnetohydro- 
dynamics (GRMHD). 

In most of these phenomena, such as for the interior dynam- 
ics of magnetized stars, or for the accretion of matter onto 
BHs, the electrical conductivity of the plasma is extremely 
high and the ideal-MHD approximation, in which the con- 
ductivity is actually assumed to be infinite, represents a very 
good approximation. In this case, the magnetic flux is con- 
served and the magnetic field is frozen in the fluid, being sim- 
ply advected with it. Following this approximation, several 
numerical codes solving the equations of general-relativistic 
ideal-MHD have been developed over the years lfH-[l3ll. By 
construction, therefore, the ideal-MHD equations neglect any 
effect of resistivity on the dynamics. In practice, however, 
even in the scenarios mentioned above, there will be spatial 
regions with very hot plasma where the electrical conductiv- 
ity is finite and the resistive effects, most notably, magnetic 
reconnection, will occur in reality. Such effects are expected 
to take place, for example, during the merger of two magne- 
tized NSs or of binary system composed by a NS and a BH, 
or near the accretion disks of AGNs, and could provide an 
important contribution to the energy losses from the system. 

The importance of resistivity effects can be easily deduced 
from the evolution of a current sheet in high but finite conduc- 
tivity. Under these conditions, several instabilities can take 
place in the plasma and release substantial amounts of energy 
via magnetic reconnection 11411 . as frequently observed, for 
example, in solar flares [15]. The study of reconnection in rel- 
ativistic phenomena is instead important to try to explain the 
origin of flares in relativistic sources, such as blazars lfl6ll and 



magnetars [17]. It is not surprising then that several groups 
have developed in the recent years numerical codes to solve 
the equations of special relativistic resistive MHD lfl8l - [24ll . 

In those scenarios that involve compact objects such as NSs 
and BHs, resistivity plays also a very important role and the 
equations of ideal MHD would not be sufficient to study them. 
In the case of general relativistic simulations of magnetized 
binary NS (BNS) and NS-BH mergers, for example, the mag- 
netic field has been up to now always confined to the interior 
of the NS, where the ideal-MHD limit is a very good approxi- 
mation l25U3lll and therefore neglecting any effect that could 
come from the magnetic field evolution in the NS magneto- 
sphere. The equations of general relativistic resistive MHD 
provide a framework that can be used to study both the re- 
gions of the domain with a high (such as the NS interior) and 
small conductivity (such as the NS magneto sphere). More- 
over, when the conductivity is set to zero, Maxwell equa- 
tions in vacuum are recovered lH9ll thus allowing for the 
study of the magnetic field evolution also well outside the NS. 
This is particularly important, since several recent works have 
claimed that the interaction of magnetic fields surrounding 
BNS and NS-BH systems may lead to strong electromagnetic 
emissions JUl, and even affect the dynamics of these systems 
(see ll33ll but also ll34ll for a different conclusion). In order 
to verify such predictions, it is therefore important to be able 
to accurately follow the dynamics of the magnetic fields in 
the region surrounding these compact binary and this cannot 
be done in the limit of ideal MHD. Last but not least, binary 
mergers are also thought to be behind the central engine of 
short gamma-ray bursts (GRBs) [[sol l35l - [37ll and the accurate 
study of the magnetic field both before and after merger could 
provide insights on current observations. 

We present the first fully general-relativistic resistive MHD 
code in a 3+1 decomposition of spacetime. We extended 
the ideal GRMHD Whisky code to include the general rel- 
ativistic version of the resistive MHD formalism presented in 
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Ref. iTToh . This new version of the Whisky code can han- 
dle different values of the conductivity going from the ideal 
MHD limit (for very high conductivities) to resistive and elec- 
trovacuum regimes (obtained respectively with low and zero 
conductivity). The code implements state-of-the-art numeri- 
cal techniques and has been tested in both fixed and dynamical 
spacetimes. In particular we show the first fully general rel- 
ativistic simulation of a magnetized NS collapse to BH using 
resistive MHD to accurately follow the dynamics of magnetic 
fields both inside and outside the NS. 

The paper is organized as follows. In SectionHIlwe describe 
the general relativistic resistive MHD equations, in Sec.HlTlthe 
main numerical methods used to solve them, and in Sec. [iVl 
our numerical tests. In Sec. [V] we summarize and conclude. 

Throughout this paper we use a spacelike signature of 
(— , +, +, +) and a system of units in which c = G = M = 
1. Greek indices are taken to run from to 3, Latin indices 
from 1 to 3 and we adopt the standard convention for the sum- 
mation over repeated indices. 



II. MATHEMATICAL SETUP 



normal to the spacelike hypersurface in a 3+1 decomposition 
of spacetime (i.e., v^n^ = 0). Notice that the time com- 
ponent is not independent due to the normalization relation 
u^u n = —1, so that 



-- W ( v l - ^ 



(5) 



where W is the Lorentz factor. 

The 3+1 decomposition of the conservation laws (0, 
© provides the evolution equations for the fluid variables 
D,U,Si, which comes from the following projections of the 
stress-energy tensor 



D 

U 
Si 

Sin 



pW, 
hW 2 



■p+-(E 2 + B 2 ), 



hW 2 v t + e ijk &B k , 
hW 2 ViVj + jijp - 



EiEj 



BiBj 



- ll0 (E 2 + B 2 ) , 



(6) 
(7) 
(8) 

(9) 



We next describe our extension of the special-relativistic 
resistive MHD formalism presented in Ref. fl9h to a general 
relativistic MHD framework. A similar (but independent) ex- 
tension has been presented recently in Il24ll . 



A. The magnetohydrodynamic equations 

The complete set of relativistic MHD equations result from 
the combination of the conservation of rest mass 

V M (^) = 0, (1) 

and the conservation energy and momentum conservation 

V V T^ V = 0. (2) 

The stress-energy tensor for a magnetized perfect fluid is 
given by 



T^v = [p{l + e)+p]u^u v +pg^ + F^F vX 



^\iv F F\ ai 



(3) 



where the rest mass density p, the specific internal energy e, 
the pressure p and the velocity describe the state of the 
fluid, and are usually referred to as the "primitive" variables. 
The pressure p is described by an equation of state (EOS) as 
a function p = p(p, e) and it is a property of the type of fluid 
considered. 

The velocity of the fluid can be decomposed as 



(4) 



where v M corresponds to the three-dimensional velocity mea- 
sured by Eulerian observers moving along a four-vector n M 



where 7^ is the usual spatial part of the metric and where 
we have introduced the specific enthalpy h = p(l + e) + p. 
The conserved rest-mass density D, the energy density U and 
the momentum Si are usually referred to as the "conserved" 
quantities since they can be shown to satisfy conservation laws 
in flat spacetimes [38] . In general, it is more convenient to 
describe the energy conservation in terms of the quantity r = 
U — D, which allows to recover the Newtonian limit of the 
energy density. 



B. The Maxwell equations 

Given a four-metric tensor g^ u , the dynamics of the electro- 
magnetic fields is described by the extended Maxwell equa- 
tions fidMi 



V V (F^ 



= —nn v 6, 



(10) 
(11) 



where F^ v is the Maxwell tensor, *F^ V is the Faraday ten- 
sor, l v is the electric current and (0, ip) are scalars to con- 
trol the constraints. In vacuum or highly magnetized plasmas, 
where the electric and magnetic susceptibilities of the medium 
vanish, the Faraday tensor can be written as the dual of the 
Maxwell tensor 



r — 2 r a(3i 



(12) 



with e^ a P 



V 



fiua/3 



I yfg and g the determinant of the four- 



metric. These tensors can be decomposed in terms of the 
electric and magnetic fields measured by an observer moving 
along a normal direction n v as: 

F^ = n»E u -n v E^ + B a np, (13) 
*F^ U = n^B v - n v B^ - e^ ap E a n^ (14) 
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Following the same decomposition, the electric current J M 
can be written as: 

I» = n^q + J M , (15) 

where q and J M are the charge density and the current for an 
observer moving along tt/\ respectively. Using these defini- 
tions and performing a 3+1 decomposition of the equations 
(fTOl) , (ITTb . (TT3T) with respect to the normal vector n^, we ar- 
rive to the following evolution equations 

(dt-CrfE* - e ijk V j (aB k ) + a~ f ij \7 j ^ = 



atrKE i -aJ\ (16) 

(d t - Cp)tj) + aV i E i =aq- (17) 
{dt-CfiB* + e ijk V j (aE k )+a 1 ij V j = 

air KB* , (18) 

(<9 t - Cp)<l> + aVi^ = -a«0, (19) 



where the scalar fields 0, -0 measure the deviation from the 
constrained solution. Their evolution equations contain damp- 
ing terms such that the constraint violations decay exponen- 
tially to zero over a timescale [Qj| HH] . 

A consequence of the Maxwell equations is the current con- 
servation 

V„ J" = 0, (20) 



which provides an evolution equation for the charge density 

(d t -Cp)q + V i (aJ i ) = aKq. (21) 

The charge density can either be computed using the evolution 
equation above or using the constraint q = ViE l . 

Finally, a relation for the current as a function of the other 
fields is needed in order to close the system. Ohm's law pro- 
vides a prescription for the spatial conduction current. We will 
consider here an isotropic scalar Ohm law 

J 1 = qv i + Wo\E l + e ijk VjB k - (v k E k )v% (22) 

where the conductivity a is chosen to be either a constant or a 
function of the rest-mass density. 



C. The full set of evolution equations 

Combining the MHD and Maxwell equations we obtain the 
following set of evolution equations, which we write in flux- 
conservative form as 





h d k {-p k ^B i + ae ik ^E j ) = 


-^B k {d k P) - a^drf, 


(23) 




V d k (-/3 k ^E i -ae ik ^B j ) = 


-^E k (d k p l ) - ay/TY^dj^ - a^jJ\ 


(24) 


d t (p - 


h d k (-p k 4> + aB k ) = -^d k p k ) 


+ B k (d k a) - ^( 7 lm d k7lm )B k - cxkcP, 


(25) 


dtip - 


V d k (-p k ^ + aE k ) = -^(d k /3 k ) 


+ E k (d k a) - ^ lm d kllm )E k + aq- a^, 


(26) 


<9i(V7<?) - 


h d k [^(-f3 k q + aJ k )} =0, 




(27) 


dt(VlD) - 


h d k [^(-p k D + av k D)] =0, 




(28) 




h d k {y^[-p k T + a(S k -v k D)}} 


= ^j(aS lm K lm - S k d k a), 


(29) 




h d k [^(-/3 k Si + aS k t )] = ^' 


\s lm d illm + S k di(3 k - (t + D)d ia . 


(30) 



r 



III. NUMERICAL SETUP 

This new version of the Whisky code implements sev- 
eral numerical methods that have been successfully used in its 
ideal-MHD version but it also implements new nu- 

merical algorithms which are instead needed in order to han- 
dle the evolution in time of the resistive MHD equations. Here 
we briefly summarize the numerical methods that are in com- 
mon with the ideal-MHD version of Whisky fTTl[27l[29il4Qll, 
while in the following section we provide a more detailed de- 
scription of the new algorithms that have been implemented. 

The evolution of the spacetime is obtained using the 



C CAT IE code, a three-dimensional finite-differencing code 
providing the solution of a conformal traceless formulation of 
the Einstein equations [40]. The general-relativistic RMHD 
equations are solved instead using high-resolution shock- 
capturing schemes (HRSC) [41]. As its ideal-MHD coun- 
terpart, also the Whisky RMHD code implements several re- 
construction methods, such as Total- Variation-Diminishing 
(TVD) methods, Essentially-Non-Oscillatory (ENO) meth- 
ods [El and the Piecewise Parabolic Method (PPM) HI. 
The Harten-Lax-van Leer-Einfeldt (HLLE) approximate Rie- 
mann solver i44ll has been used to compute the fluxes in all 
the results presented here. Since the code is based on the 
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Cactus 14511 computational framework, it can also use adap- 
tive mesh refinement (AMR) via the Carpet driver l46ll . 



A. IMEX Runge-Kutta Methods 

The general-relativistic RMHD equations in high- 
conductivity media contain stiff terms which make the time 
evolution with an explicit time integrator very inefficient, if 
not impossible. The prototype of the stiff system of partial 
differential equations can be written as 



d t U = F(U) + Jii(U), 



(31) 



where e = l/a > is the relaxation time. In the limit of 
e — » oo, the second term on the right-hand side of Eq. (PIT) 
vanishes and the system is then hyperbolic with a spectral ra- 
dius Ch (i.e., with Ch being the absolute value of the maxi- 
mum eigenvalue). In the opposite limit of e — >• the first 
term on the right-hand side of Eq. (|3TT) vanishes and the sys- 
tem is clearly stiff, since the timescale e of the relaxation (or 
stiff term) i?(U) is very different from the speeds Ch of the 
hyperbolic (or non-stiff) part F(XJ). 

Stiff systems of this type can be solved efficiently by a com- 
bination of implicit and explicit time integrators. In particular, 
the IMEX Runge-Kutta scheme consists in applying an im- 
plicit discretization to the stiff terms and an explicit one to the 
non-stiff terms. When applied to the system (|3TT) it takes the 
form El 

i-l 



u (i) =U n + A^a^-F(U^), 
N 1 

3 = 1 £ 
N N 

At WiF(U® ) + At^ oUi-R(U {i) ) , 



u n+1 = u n 



where are the auxiliary intermediate values of the Runge- 
Kutta time integrator. The matrices A = (dij), dij = for 
j > i and A = (a^), are N x N matrices such that the 
resulting scheme is explicit in F and implicit in R. An IMEX 
Runge-Kutta scheme is characterized by these two matrices 
and the coefficient vectors Cbi and c^. Since the simplicity 
and efficiency of solving the implicit part at each step is of 
great importance, it is natural to consider diagonally-implicit 
Runge-Kutta (DIRK) schemes for the stiff terms, i.e., (a^ = 
for j > i). The matrices of coefficients are reported in Table 

m 

Our approach to the solution of the potentially stiff set 
of general-relativistic RMHD equation consists therefore in 
the use of the IMEX RK method introduced above. For 
the particular set of equations (f23l)-(f30lh the evolved fields 
can be split into stiff terms V = {E 1 } and non-stiff terms 
W = {B*,^,0,g,r,5 < ,D}. 

The evolution of the electric field (f24l) can become stiff de- 
pending on the value of the conductivity a = 1 je in the Ohm 



TABLE I: Tableau for the explicit (left) implicit (right) IMEX- 
SSP3(4,3,3) L-stable scheme 
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1/2 
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c 1/2 


-b-c- 


a a 







1/6 


1/6 


2/3 



a = 0.24169426078821 , b = 
c = 0.12915286960590 



0.06042356519705 , 



law (l22l) . Its right-hand- side can be splitted in potentially- stiff 
terms and regular ones 



d t {yfiE % ) = F % E + Re, (33) 
where we have introduced the factor 1 je on the definition of 



R l E and 



F E = -d k [-f3 k SYE*-ae ik ^Bj]-^E k (d k F)- 



R\ 



(v k E k y] . 



(34) 
(35) 



In order to evolve this system numerically, the 
fluxes {F T ,F S i,Fr)} have to be computed at each 
timestep. This implies that the primitive quantities 
{p, p, v l , E % , B 1 } have to be recovered from the con- 
served fields {D, r, S\ E\ ^ B 1 }. With the 
exception of very simple EOSs, this recovery cannot be 
done analytically and it is instead necessary to solve a set of 
algebraic equations via some root-finding iterative procedure, 
which we will describe below. 

Before that, we note that the solution of the conserved quan- 
tities {D, r, S\ B 1 } at timet = (n+1) At is obtained by 
simply evolving the equations (l28k (l30b , d23l) . However, the 
same procedure for the electric field leads only to an approx- 
imate solution {E 1 } containing only the explicit terms. The 
full solution, involving also the potentially stiff terms, requires 
therefore the inversion the implicit equation (l24k which de- 
pends on the velocity v % and the fields {B\ E 1 }. In the case 
of the scalar Ohm law (|22l) . the stiff part is linear in E l , so a 
simple analytic inversion can be performed 

E i = M~ 1 (v j ) [E* + a S E (v\B j )}, (36) 

where a = an At a W a and the inversion matrix is given by 



M= 



1 + a(l — v x v x ) —a(v y v x ) 



-a(v x v y ) l + a(l-v y v y ) 
-a(v x v z ) 



-cr(v z v x ) 
-a(v z v y ) 



cr(v y v z ) 1 + cr(l - v z v z ) 

(37) 



The recovery procedure is similar to the one presented in 
Ref. fl9ll and can be summarized in the following steps: 
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1. Consider an initial guess for the electric field. Some 
possible options are: its value in the previous timestep, 
its approximate value in the current time step E\ or the 
ideal MHD value E % = —e^ k VjBk, where Vj is the 
velocity in the previous time level. 

2. Subtract the electromagnetic field contributions from 
the conserved fields, namely compute 



f = T -\{E* + B% 
Si = Si CijkE^ B . 



(38) 
(39) 



3. Perform the recovery as in the non-magnetized case: 
The EOS can be used to write the pressure as a function 
of the conserved quantities and the unknown x = hW 2 , 
so that the definition of r can be written as 



i 



w 2 r 
r-r„ 



i 



r(r P 



D K lw 



YW 



12? 



(40) 



which must vanish for the physical solutions. Here T p 
and T are the adiabatic indices corresponding to an ideal 
gas and a poly tropic EOS, respectively, while K is the 
polytropic constant. By setting T = 1 we recover the 
simple polytropic EOS, while the ideal EOS can be re- 
covered by setting Y p = T. 

4. A solution of the function f(x) = can be found 
numerically by means of an iterative Newton-Raphson 
solver. The initial guess for the unknown x is given by 
the previous time step. 

5. After each step of the Newton-Raphson solver, update 
the values of the fluid primitives 



Vi = ^, W 2 = 



S 2 



r- i / x 



v 



{w 2 ~ p ) 



D 



r \w 2 r ) 1 r(r p -i) " 

and then invert the electric field according to (l36l) . 



(41) 



(42) 



IV. NUMERICAL TESTS AND RESULTS 

In this extended Section we report the numerical results ob- 
tained in one-, two- and three-dimensional tests, which con- 
firm that our implementation is correct and provides the ex- 
pected results in a large range of conductivities. More specif- 
ically, the one-dimensional tests involve: i) a large- amplitude 
circularly-polarized (CP) Alfven wave to validate that our im- 
plementation matches the ideal-MHD results in the high con- 
ductivity regime; ii) the evolution of a self- similar current 
sheet, which tests our implementation in the intermediate con- 
ductivity regime; Hi) a collection of shock- tube tests involving 
a range of uniform and non-uniform conductivities. In these 
particular tests we also examine the zero-conductivity regime, 
where the electromagnetic fields are expected to follow the 
vacuum Maxwell equations and hence behave as propagating 
waves. 

Following the one-dimensional tests, we then present two 
and three-dimensional tests, which include the standard cylin- 
drical and spherical explosion tests, which we consider in the 
case of very large conductivities in order to test the ideal- 
MHD limit of our equations. Finally, we have performed three 
different sets of simulations involving spherical magnetized 
stars in general relativity. The first setup consists in a spherical 
(TOV) star with prescribed magnetic fields confined initially 
in the interior of the star. The second set involves the evo- 
lution of a magnetized star with initial data generated by the 
LORENE library and having a dipolar magnetic field that ex- 
tends also outside the star. As a conclusive three-dimensional 
test we consider the gravitational collapse of a nonrotating star 
to a black hole, where the initial data is again generated by the 
LORENE library H]. 

With the exception of the collapsing star, where we have 
used a polytropic EOS, all simulations reported here have em- 
ployed an ideal gas (T-law) EOS 



v = pe(r - 1) , 



(43) 



with T = 2 for the one-dimensional tests and T = 4/3 for the 
two and three-dimensional tests. In addition, for the evolu- 
tion of the stable magnetized stars we have adopted a T = 2. 
As mentioned above, the collapse of the unstable magnetized 
star has been followed using a polytropic EOS, p = Kp F , 
with T = 2. Finally, to ensure a divergence-free magnetic 
field with the our hyperbolic divergence-cleaning approach, 
we have set the damping coefficient n to be one everywhere. 



6. Iterate the steps 2.-5. until the difference between two 
successive values of x and the electric field fall below a 
given threshold, usually of the order of 10 ~ 10 . 

This procedure converges quickly in the high-conductivity 
regions if the ideal MHD solution is chosen as an initial guess, 
and in the intermediate conductivity regions if the initial guess 
is given by the approximate electric field E % . In general, < 5 
iterations are sufficient for intermediate conductivities, while 
< 70 iterations are usually necessary in the regions with high 
conductivity. 



A. One-dimensional Test Problems 

1. Circularly Polarized Alfven waves 

The present test has been discussed in detail in Ref . II 1 (Jl and 
it computes the propagation of a large- amplitude circularly- 
polarized Alfven wave through a uniform background mag- 
netic field Bq. For the purpose of this test, we set a very high 
conductivity a = 10 6 in order to recover the ideal-MHD limit. 
Since the propagating wave is expected to be the advected ini- 
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FIG. 1: Circularly-Polarized Alfven wave. B y component 
of the magnetic field for three different resolutions Ax = 
{1/50, 1/100, 1 /200}, together with the exact initial solution (black 
solid line). Clearly, the numerical solution provided by the resistive 
MHD implementation and the exact one overlap for a uniform con- 
ductivity a — 10 6 and the highest resolution. 



FIG. 2: Self-similar current sheet. B y component of the magnetic 
field at the initial t = 1 and final time t — 10. The exact solution at 
t = 1 is shown with a dashed blue line. The solution given by the 
analytic expression d4Tb at t = 10 (black solid line) is indistinguish- 
able from the numerical solution obtained form the resistive MHD 
equations (red dashed line). 



tial profile, it is convenient to apply periodic boundary con- 
ditions and compare the evolved profile after one full period 
with the initial one, in order to check the accuracy of our im- 
plementation. 

In particular, we consider a CP Alfven wave with a normal- 
ized amplitude t]a traveling along positive x-axis, in a uni- 
form background magnetic field Bo with components 

B % = {B , tjaBo cos[k(x - VAt)],rjABo sm[k(x - VAt)]}. 

(44) 

For simplicity, we take v x = and write the remaining veloc- 
ity components as 



-v A By/B , 



-v A B z /B , 



(45) 



where 



2Bl 




KP h + B$(l+ V \) 

(46) 

By setting p = p = tja = 1 and Bq = 1.1547, we fix the 
Alfven velocity to va = 0.5. Therefore, in a computational 
domain centered at x = with x G [—0.5, 0.5], we expect 
the wave to return to its initial position after one full period 
t = L/va = 2. The comparison of the numerical solution 
with the initial condition (l44b at t = gives us a measure of 
the error. 

In principle, the resistive MHD formalism would allow us 
to recover the ideal-MHD limit only for an infinite conduc- 
tivity. In practice, however, the use of a conductivity as large 
as a = 10 6 is sufficient to obtain a solution that converges 
to the ideal-MHD one with increasing resolution. As a re- 
sult, we have chosen to perform simulations with a uniform 



conductivity of a = 10 6 , using the following resolutions: 
Ax = {1/50,1/100,1/200}. 

In Fig. [T] we show the component B y at time t = 2, cor- 
responding to one full period. By superimposing the results 
at t = 2 with the initial data at t = 0, it is evident that the 
numerical solution of the resistive MHD equations tends to 
the ideal-MHD exact solution for a high-enough conductiv- 
ity and resolution. We have used both a linear reconstruction 
with monotonized-central (MC) slope limiter and the second 
order PPM reconstruction. The numerical solution converges 
to the exact one at second order when using PPM reconstruc- 
tion and at second order with the linear reconstruction, exactly 
the same convergence rates than with the original ideal MHD 
system implemented in WhiskyMHD. 



2. Self-similar Current Sheet 

We next considera a test that involves the evolution of a self- 
similar current sheet, as proposed in Ref. lfl8ll . This setup is 
useful for testing codes which solve the resistive MHD equa- 
tions with a moderate conductivity regime, which we set to be 
a = 100. 

In practice, the initial data consists in a magnetic field 
solely in the ^-direction which changes sign in a thin cur- 
rent layer. Provided that the initial solution is in equilibrium 
(i.e., the pressure and density are constant, and the velocity 
is zero) and that the magnetic pressure is much smaller than 
the fluid pressure everywhere, then the evolution of the mag- 
netic field is given by the simple diffusion equation dfB y — 
(1/cr) d%B y = 0, which will be responsible for the diffu- 
sive expansion of the current layer in response to the physical 
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FIG. 3: Shock-tube Tests. B y component of the magnetic field at 
t = 0.4 for different resolutions Ax = {1/100, 1/200, 1/400}. 
The highest resolution Ax = 1/ 400 matches the exact ideal-MHD 
solution remarkably well. 

resistivity (we are also assuming that E % = = d t E % ). Un- 
der these simplified assumptions, the analytical solution of the 
diffusion equation is given, for t > 0, by 

B"(M)=B Erf Q^l) > ( 47 ) 

where £ = t/x 2 and Erf is the error function. Clearly, as the 
evolution proceeds, the current layer expands in a self- similar 
fashion. 

Following OjlQjl], we use as initial data the analytic so- 
lution (l47l) at t = 1 and set the density and pressure to be 
uniform with p = 1 and p = 50 respectively, while keep- 
ing the components of the electric field and velocity to zero 
initially. In our calculations we have used a computational 
domain with extents x = y = z G [— 5, 5] with a resolution 
of Ax = 1/200. Furthermore, a linear reconstruction method 
was adopted with the further application of the MC limiter. 

In Fig. [2] we present the results we obtained by solving 
numerically the resistive MHD equations and the compari- 
son with the exact solution (1471) at t = 10 (black solid line). 
Clearly, the numerical solution (red dashed line) is indistin- 
guishable from the analytic one, thus providing convincing 
evidence that the code can accurately describe resistive evolu- 
tions with intermediate values of the conductivity. 



3. Shock-Tube Tests 

We next consider the numerical solution of the standard of 
Brio and Wu shock-tube test [49] as adapted for its MHD im- 
plementation and using either a variety of uniform or space- 
dependent conductivities parameterized by the reference con- 
ductivity ctq. More specifically, the left (L) and right (R) 



FIG. 4: Shock-tube Tests. B y component of the magnetic field for 
conductivities a = {0, 10, 10 2 , 10 3 , 10 6 } at t = 0.4 and resolution 
Ax = 1/200. For ao = the magnetic field is governed by a wave- 
like equation, corresponding to the solution of the Maxwell equations 
in vacuum. 

states are initiall y se parated by a discontinuity at x = 0.5 
and are given by 113011 

(p L , PL, B\) = (1.0, 1.0, 0.5), 
(pr,Pr, B y R ) = (0.125, 0.1, -0.5), 

while all other variables are set to zero. The ideal-MHD evo- 
lution of the aforementioned setup with B x = leads to 
two fast waves, one rarefaction propagating to the left and a 
shock propagating to the right of the discontinuity. The solu- 
tion of this test in the ideal-MHD limit exists and is found in 
the exact ideal-MHD Riemann solver provided by Ref. [50]. 
For the rest of the one-dimensional tests, any comparison be- 
tween the solution of the resistive MHD equations in the high- 
conductivity regime and the exact solution of the ideal-MHD 
equations is performed with data obtained from the publicly 
available code [50]. All tests have been performed employing 
a linear reconstruction method with further application of the 
MC slope limiter. 

As a first setup of our shock- tube tests, we consider the 
case of a uniform high conductivity a = ao = 10 6 and, in 
analogy with the Alfven-wave test in the high-conductivity 
regime, we verify that the solution of the coupled Maxwell- 
Hydrodynamics equations tends to the ideal-MHD exact so- 
lution Il50ll as the resolution is increased. Figure [3] reports 
the magnetic field component B y at t = 0.4 for the three 
resolutions Ax = {1/100, 1/200, 1/400} considered. The 
high-resolution result matches the exact ideal-MHD solution 
so well that is difficult to distinguish them, thus providing the 
first evidence that our implementation is robust also in the 
presence of discontinuities. 

As a second setup of the shock- tube tests, we consider the 
case in which the conductivity is still uniform in space, but 
of different strength. In particular, we perform the same test 
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FIG. 5: Shock-tube Tests. The left panel shows the conductivity profile at £ = 0.4 for non-uniform conductivity with different power laws, 
i.e., 7 = {0,6,9,12}. The 7 = case corresponds to the high-conductivity regime of the resistive MHD equations. The right panel reports 
instead the B y component of magnetic field for the same initial conditions as in the left one. The leftmost region tends to the ideal MHD 
solution, while the rightmost tends to the vacuum solution for 7=12. 



for a = {0, 10, 10 2 , 10 3 , 10 6 }, while keeping the resolution 
fixed at Ax = 1/200. Figure [4] reports different solutions of 
the magnetic-field component B y given by the resistive MHD 
equations with different values of cfq . It is important to note 
here that the solutions change smoothly from the ideal-MHD 
solution computed for cfq = 10 6 , to the wave-like solution 
for do = 0, which corresponds to the propagation of a dis- 
continuity at the speed of light, corresponding to a solution 
of the vacuum Maxwell equations. The ability of treating the 
two extreme behaviours of the Maxwell-MHD equations via 
a resisitive treatment is an essential feature of our approach 
and a fundamental one in the description of the dynamics of 
magnetized binary neutron stars. 

As a final setup our of our suite of shock-tube test, we have 
considered the same initial data but now prescribed a non- 
uniform conductivity given by the expression 

cj = a D^ , (48) 

where 7 is an integer exponent we vary in the range 7 G 
[0,12]. Thes prescription above introduces nonlinearities with 
respect to the conserved rest-mass density D and provides 
an intuitive way of tracking the dense fluid regions. It leads 
to low values of the conductivity in places were the plasma 
is tenuous and high values in more dense regions, which 
will prove very useful later on when evolving magnetized 
stars. However, this prescription is far from being realis- 
tic and normally a more general conductivity prescription 
a = a(D,T,E) is to be seeked starting from micro-physical 
considerations. 

Following [19], we adopt the same initial data as before, 
however this time we change the exponent 7 of Eq. ([48b while 
maintaining the value of conductivity to ao = 10 6 . 

The results of this last test are reported in the left panel of 
Fig- 13 which show the profile of the conductivity at £ = 0.4 



for different values of the power-law exponent, i.e., 7 = 
{0, 6, 9, 12}. Clearly, the conductivity follows the evolution 
of the rest-mass density, with a left-going rarefaction wave 
and right-going shock. It is interesting to note that our ap- 
proach is able to track even very large variations in the con- 
ductivity, with jumps as large as eleven orders of magnitude 
across the computational domain. The right panel of Fig.[5j on 
the other hand, reports instead the magnetic field-component 
B y at £ = 0.4 for the same initial conditions. As imposed by 
Eq. d48b , the solution in the leftmost part of the computational 
domain, where the rest-mass density is very high, is controlled 
by a very high conductivity, which tends to <Jo = 10 6 . In turn, 
this implies that the solution for the magnetic field should ap- 
proach the ideal-MHD limit in that region. On the other hand, 
in the rightmost region, where the rest-mass density is very 
low, the conductivity is correspondigly small and tending to 
zero for high values of 7. In such regions, therefore, the mag- 
netic field is expected to behave as a wave, thus explaining the 
appearance of a moving peak for 7 = 12. 

Overall, this suite of shock-tube tests, demonstrates that 
our numerical implementation is able to treat accurately both 
uniform and non-uniform conductivity profiles in one dimen- 
sional tests, independently of the steepness of the profiles and 
even in the presence of shocks. 

B. Multidimensional Tests 

We now focus on multidimensional tests that involve 
shocks in several directions, such as the two-dimensional 
cylindrical explosion and the three-dimensional spherical ex- 
plosion test suggested in Ref.[l]. Despite the fact that there 
is no analytical solution for any of these tests, even in the 
ideal-MHD case, the symmetries of the problem can be of 
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FIG. 7: Left Panel: One-dimensional cuts along the ^-direction and at t = 4.0 of the the pressure. The black dashed line corresponds to the 
resistive code (the WhiskyRMHD code), while the blue dotted line corresponds to the ideal-MHD code, (the WhiskyMHD code). Right Panel: 
The same as in the left panel but for the Lorentz factor. 



great help in verifying that the numerical implementation is 
correct and that it preserves the expected symmetries. Our ap- 
proach in these tests will be therefore that of comparing the 
solution of the same multidimensional test as obtained with 
the ideal-MHD code presented in fTTIl and our new resistive 
WhiskyRMHD code in the limit of very high conductivities. 
The initial electric field is computed in such a way that it sat- 
isfies the ideal-MHD condition, i.e., E % = —e^ k VjBk, and all 
the tests have been performed adopting a linear reconstruction 
method and the minmod slope limiter. 



1. Cylindrical Blast-Wave 

In the two-dimensional cylindrical blast- wave problem, we 
adopt a square domain with 200 grid cells per direction, in a 
range of (—6.0, 6.0) x (—6.0, 6.0). The setup of the problem 
consists of three regions. The innermost region with < r < 
0.8, for which the pressure and the density are set to p = 1, 
p = 0.01 respectively, the intermediate region which extends 
from 0.8 < r < 1.0 where r = (x 2 + y 2 ) 1 ^ 2 both the pressure 
and the density exponentially decrease, and the outermost re- 
gion which is filled with an ambient plasma with p = 0.001, 
p = 0.001 and occupies the domain 1.0 < r < 6.0. The 
initial magnetic field is along the ^-direction with an initial 
magnetic field strength of Bq =0.05. 
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The numerical results are presented in Fig. [6j where we 
show that the magnetic field solution is regular everywhere 
and that there are no visible artifacts that could indicate a 
possible symmetry error in our implementation. Furthermore, 
when one-dimensional cuts of the resistive solution are plot- 
ted against the ideal-MHD solution obtained with the code 
presented in ifTTIl . the agreement is extremely good (this is not 
shown in Fig. [6]). 



2. Spherical Blast-Wave 

In the three-dimensional spherical blast- wave problem, the 
grid structure is similar, but the domain is now within the 
ranges (-6.0,6.0) x (-6.0,6.0) x (-6.0,6.0). The prob- 
lem setup consists of the same three regions as in the cylin- 
drical blast wave problem, although here the radius r refers to 
the spherical-polar radial coordinate, and not to the cylindrical 
radius, i.e., r = (x 2 + y 2 + z 2 ) 1 / 2 . 

The corresponding solution of the spherical blast-wave 
problem in the (x, y) plane is essentially identical to the one 
already reported in Fig. [6] and for this reason we do not show 
it here. What we do show in Fig. [7] however, are one- 
dimensional cuts along the z-direction of the pressure p and of 
the Lorentz factor W as computed with the ideal-MHD code 
(blue dotted line) and the resistive MHD code (black dashed 
line). This comparison, which is not expected to be exact 
given that the resistivity is large but not infinite, provides con- 
vincing evidence of the ability of our implementation to accu- 
rately describe higher-dimensional discontinuous flows in the 
high-conductivity regime. 



C. Nonrotating Magnetized Stars 

In the following Section we present the numerical results 
obtained from the evolution of nonrotating spherical stars in 
the presence of electromagnetic fields and for a variety of con- 
ductivities. In order to accurately model both the interior and 
the exterior of the star, we prescribe a spatial dependence of 
the electrical conductivity such that the ideal-MHD limit is 
recovered in the deep interior of the star (which is expected 
to be an excellent conductor) and such that the electro vacuum 
limit is recovered outside the star, where the density and the 
isotropic conductivity is expected to be negligibly small. 

This behaviour can be easily achieved assuming that the 
conductivity tracks the (conserved) rest-mass density, thus in- 
suring a smooth transition between the two regimes. In prac- 
tice, we have experimented with functional prescriptions of 
the type 

a = cr max [(1 - D 

atmo 

/D),0] 2 , (49) 

where a ~ ao is the conductivity in the regions of large rest- 
mass density (a — ao at the stellar center) and a = in 
the atmosphere, where we set the conserved rest-mass den- 
sity to its uniform value D = D atmo . In our calculations we 
normally set do = 10 6 and D a t mo to be about ten orders of 



magnitude smaller than the value of D at the center of the 
star. Furthermore, in the atmosphere we set the fluid veloc- 
ity to zero and since a = there, the electric and magnetic 
fields are evolved via the Maxwell equations with zero cur- 
rents (electro vacuum). 

This non-uniform conductivity prescription allows us to 
provide effective boundary conditions at the surface of the 
star for the exterior electrovacuum solution similar to those 
in Refs. [EliEl], but without the limitations of using an ana- 
lytical solution for the interior of the star or the further com- 
plications of finding a suitable matching between the electro- 
magnetic fields of the interior ideal-MHD solution and the ex- 
terior one. All the simulations reported hereafter have been 
performed adopting the PPM reconstruction scheme, for rel- 
ativistic stars whose initial properties are summarized in Ta- 
bleUD 

1. Stable Star with confined magnetic fields 

For the sake of simplicity, we consider as initial data spher- 
ical stars in equilibrium to which a poloidal magnetic field 
confined to the stellar interior is superimposed (see, e.g., ll53l - 
D5I1 ). While the hydrodynamical quantities are consistent so- 
lutions of the Einstein equations, the magnetic field is added 
a-posteriori, with a consequent violation of the constraint at 
the initial time. In practice, however, this violation is always 
very small, even for the largest fields, and is quickly domi- 
nated by the violations introduced by the standard evolution. 

The toroidal vector potential that generates the poloidal in- 
terior magnetic field is expressed as fllll 

= r 2 max [A b (P - P cut ), 0] 2 , (50) 

where P cut is about 4% the central pressure P c . The the star, 
initially computed with a poly tropic EOS with T = 2, K = 
100, has a gravitational mass M = 1.4OM and is endowed 
with a poloidal magnetic field of strength B c = 10 12 G at the 
center of the star and f3 = p m &g/p = 4.49 x 10 -13 , with p mag 
the magnetic pressure. The magnetic field in the atmosphere 
is initially zero. For all of the evolutions presented hereafter 
we have used an ideal-fluid EOS with T = 2. 

We first examine the evolution of the magnetized star 
in the fixed spacetime of the initial solution (Cowling- 
approximation). In the left panel of Fig. [8] we show with thin 
solid, dashed and dotted lines the evolution of the central rest- 
mass density normalized to its initial value p c ,o in thin colored 
lines. The tests were performed using three spatial resolutions 
of Ax = {0.4, 0.3, 0.2} km, corresponding respectively to 
TV = {80, 120, 160} points across the finest AMR grid, which 
extends up to R out = ±18 km. As customary in this type 
of tests, stellar oscillations are triggered by the truncation er- 
ror and their amplitude decreases as the numerical resolution 
is increased. The importance of the test lays therefore in the 
calculation of the eigenfrequencies of the oscillations, which 
we find to be in very good agreement with those computed 
via perturbative analyses (not shown here) and with other hy- 
drodynamics and ideal-MHD codes fTTll56ll . In addition, a 
comparison with the ideal-MHD code 1 1 1] shows a very good 
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Star type 


Madm [M©] 


M b [M Q ] 


R eci [km] 


K 


r 


B c [G] 


# levels 


N 


-/Vstar 


Rout [km] 


Confined fields 


1.40 


1.51 


12.00 


100.0 


2 


10 12 


4 


80, 120, 160 


56,80,112 


142 


Extended fields 


1.33 


1.37 


32.56 


372.0 


2 


2.4 x 10 14 


4 


120 


84 


355 


Unstable model 


2.75 


2.89 


16.30 


364.7 


2 


5 x 10 15 


5 


272 


216 


241 



TABLE II: Properties of the magnetized star models used in the simulations. The columns report: the ADM and baryon masses in units of 
solar masses Madm and Mb respectively, the circumferential equatorial radius of the star in kilometers i? eq , the polytropic constant K, the 
polytropic index T, the value of the magnetic field in Gauss at the center of the star B c , the number of refinement levels, the number of 
gridpoints on the finest level N, the number of gridpoints across the star N star for the different resolutions considered, the computational grid 
outer boundary in kilometers R out . 




FIG. 8: Left Panel: Evolution of the central rest-mass density of a nonrotating magnetized star for both the Cowling approximation (C; thin 
lines) and a dynamical spacetime (D; thick lines). Different line types mark different resolutions: dashed light blue Ax = 0.4 km, dotted dark 
blue Ax = 0.3 km, continuous black Ax = 0.2 km. Middle Panel: The same as the left one but for the central magnetic field. Right Panel: 
The same as the middle one but different values of the conductivity ctq. All lines refer to a resolution of Ax = 0.2 km. 



agreement in the evolution of the rest mass density, indicat- 
ing that the oscillations are tracked correctly by our resistive 
MHD implementation. 

We next examine the same scenario, but in a fully dynam- 
ical spacetime and find also in this case a very good agree- 
ment with the ideal MHD solution. Still in the left panel of 
FigIS]we report with thick solid, dashed and dotted lines the 
evolution of the central rest-mass central density in a dynam- 
ical spacetime for different resolutions. As well known from 
perturbation theory, the eigenfrequencies of oscillations are in 
this case lower but what is relevant to note is that the secular 
evolution in both the fixed and dynamical spacetimes are very 
similar, with variations in the central density that is less than 
a couple of percent over tens of dynamical timescales. 

The middle panel of Fig. [8] displays instead the evolution 
of the central value of the magnetic field, where lines of dif- 
ferent color refer to different resolutions, while the thickness 
marks whether we are considering a fixed or a dynamical 
spacetime (thin for the Cowling approximation and thick for a 
full general-relativistic evolution). The corresponding power 
spectral density is shown in Fig. [9j where different line types 
refer to different resolutions and the dotted vertical lines mark 
the eigenfrequencies obtained from linear perturbation theory. 
The match between the numerical and perturbative results is 
clearly excellent and the differences in the fundamental mode 
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FIG. 9: Power spectral density of a full general-relativistic evolution 
of the central rest-mass density for a stable star with confined mag- 
netic fields. Different line types refer to different resolutions. Shown 
with dotted vertical lines are the eigenfrequencies obtained from lin- 
ear perturbation theory. 
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FIG. 10: Two-dimensional cuts on the (x, z) plane of the solution the rest-mass density (colorcode from white to red) and of the magnetic 
field lines and at times t = 0, 9.88, and 18.59 ms. The evolution refers to a nonrotating star in a dynamical spacetime. Note that although the 
magnetic field is contained in the star initially, it diffuses out as a result of numerical and physical resistivity. 



at the highest resolution are < 0.5%. 

We note that, as for the central rest-mass density, the evolu- 
tion of the central magnetic field is accompanied by a secular 
drift towards lower values, and this is simply the result of the 
intrisic numerical resistivity (we recall that these tests have 
been performed with the resistive code but for very large con- 
ductivities and hence in a virtual ideal-MHD regime). Clearly, 
the numerical resistivity decreases with resolution and this is 
exactly what the behaviour in the middle panel shows. It is 
interesting to note that while with sufficient resolution the re- 
sistive losses saturate to about 20% of the original magnetic 
field over ~ 12 ms, these can be very large for low resolution 
and dissipate up to ~ 85% of the initial magnetic field over 
the same time-span. These numerical resistive losses should 
be compared with the ones introduced instead by the physi- 
cal resistivity and which can of course be much larger. This 
is shown in the right panel of Fig. [8] which is the same as 
the middle one, but where we have used the highest resolution 
(i.e., Ax = 0.2 km) and varied the strength of the physical re- 
sistivity from do = 10 6 to ao = 10 2 . Because the fluid veloc- 
ities are essentially zero at this resolution, the magnetic-field 
evolution follows a simple diffusion equation with a Ohmic 
decay timescale which scales linearly with 1/cr. This is in- 
deed what shown in the right panel of Fig. [8] where, after the 
initial transient, the solution settles to an exponential decay 
and where the magnetic field can be reduced of almost two 
orders of magnitude over 12 ms in the case of ao = 10 2 . 

Finally, we show in Fig. [TO] two-dimensional cuts on the 
(x, z) plane of the rest-mass density (shown in a colorcode 
from white to red) and of the magnetic field lines for an os- 
cillating star; the three panels refer to times t = 0, 9.88, and 
18.59 ms, respectively. It is important to remark that although 
we start with a magnetic field that is initially confined inside 
the star, the inevitable presence of a small but finite numerical 
resistivity and our choice of a nonzero physical conductivity 
near the surface of the star [we recall that our conductivity 
follows the profile given in Eq. (09])], induce a slow but con- 
tinuous "leakage" of the magnetic field, which leaves the star 
and fills the computational domain. Because the external mag- 
netic field is essentially with a zero divergence and with a van- 



ishingly small Laplacian (we recall that in the stellar exterior 
the resistivity is zero and the Maxwell equations tend to the 
those in vacuum), it is to a very good approximation a poten- 
tial field, as shown by the clean dipolar-like structure. Clearly, 
the Ohmic diffusion timescale increases with resolution and 
therefore the relaxation of the magnetic field to a stationary 
dipolar-like structure takes place on longer timescales for the 
high-resolution simulation. 



2. Stable Star with extended magnetic fields 

We next consider initial data for a spherical magnetized star 
with a poloidal magnetic field extending outside the star, as 
generated by the Mags tar code from LORENE library l48ll . 
The external magnetic field is dipolar and is computed by 
solving the Maxwell equations in vacuum, with boundary con- 
ditions given by the interior poloidal magnetic field. This so- 
lution is fully consistent with the Einstein equations and it 
provides accurate measurements of the stellar deformations in 
response to either rapid rotation or large magnetic fields ll57ll . 
More specifically, we have considered a nonrotating star mod- 
elled initially as poly trope with T = 2 and K = 372, hav- 
ing a gravitational mass M = 1.33 M©, and endowed with 
a poloidal magnetic field of strength £? c = 5xl0 15 G. The 
magnetic field in the atmosphere is given by the electrovac- 
uum solution, which has a dipolar structure. The evolutions 
have been carried out in a computational domain with outer 
boundary at R out = 355 km and a resolution of Ax = 0.7 
km, corresponding to 60 points covering the positive part of 
finest grid which extends up to 44 km. 

Figure [TT] displays in its left and middle panels two- 
dimensional cuts on the (x, z) plane of the rest-mass density 
(shown in a colorcode from white to red) at the initial and fi- 
nal times, i.e., t = ms and t = 37 ms. A rapid comparison 
among the three panels clearly shows the ability of the code 
to reproduce stably the evolution of this oscillating star also 
when the magnetic field extends in its exterior. The right panel 
of Figure[TTl on the other hand, shows in its top part shows the 
evolution of the magnetic flux computed across a hemispheric 
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FIG. 11: Left and Middle Panels: Evolution of the magnetic field lines displayed at times t = Oms and t = 37 ms. The rest mass density is 
also shown with purple-red-yellow colors. Right Panel: The top part shows the evolution of the magnetic flux computed across a hemispheric 
surface at a radius r = 135 km, while the bottom part shows the power spectral density of the rest-mass density (black solid line) and of the 
magnetic flux (blue dotted line). 



surface at a radius r = 135 km, which also shows signs of 
oscillations. We have computed the power spectrum of these 
oscillations and compared it with the corresponding one ob- 
tained for the central rest-mass density. The results of this 
comparison are shown in the bottom part of the right panel, 
with a black solid line referring to the rest-mass density a blue 
dotted line to the magnetic flux. The very good agreement 
between the two implies that the oscillations observed in the 
magnetic flux are essentially triggered by the oscillations in 
the rest-mass density. 



3. Magnetized Collapse to a Black Hole 

Our final and most comprehensive test is represented by 
the collapse to a BH of a magnetized nonrotating star. This 
is more than a purely numerical test as it simulates a pro- 
cess that is expected to take place in astrophysically realis- 
tic conditions, such as those accompanyin g th e merger of a 
binary system of magnetized neutron stars |2a,[27|], or of an 
accreting magnetized neutron star. The interest in this pro- 
cess lays in that the collapse will not only be a strong source 
of gravitational waves, but also of electromagnetic radiation, 
that could be potentially detectable (either directly or as pro- 
cessed signal). The magnetized plasma and electromagnetic 
fields that surround the star, in fact, will react dynamically 
to the rapidly changing and strong gravitational fields of the 
collapsing star and respond by emitting electromagnetic radi- 
ation. Of course, no gravitational waves can be emitted in the 
case considered here of a nonrotating star, but we can never- 
theless explore with unprecedented accuracy the electromag- 
netic emission and assess, in particular, the efficiency of the 
process and thus estimate how much of the available bind- 
ing energy is actually radiated in electromagnetic waves. Our 
setup also allows us to investigate the dynamics of the elec- 
tromagnetic fields once a BH is formed and hence to assess 
the validity of the no-hair theorem, which predicts the expo- 
nential decay of any electromagnetic field in terms of Quasi 
Normal Mode (QNM) emission from the BH. 



Ours is not the first detailed investigation of this process 
and relevant previous studies are that in Ref. Dill and the 
more recent one in Ref. i52ll . However, our approach differs 
from previous ones in that it correctly describes the gravita- 
tional dynamics of a collapsing fluid (the semianalytical work 
in Ref. Dill , in fact, considered the more rapid collapse of a 
dust sphere, for which the Oppenheimer- Snyder (OS) analytic 
solution can be used ll58ll ) and does not require any match- 
ing of the solution near the stellar surface (the fully relativis- 
tic work in Ref. ll52ll had to resort to an ingeniuous match- 
ing between the interior ideal-MHD solution and a force-free 
one in the magneto sphere), leaving the complete evolution of 
the electromagnetic fields to our prescription d49b of a non- 
uniform conductivity. Indeed, our solution is expected to be 
exactly the same as the force-free one except in regions where 
B 2 — E 2 < and an anomalous resistivity appears. Since 
we can handle accurately such resistive regions, this test il- 
lustrates the capabilities of our resistive implementation and 
serves as a more realistic approach to this astrophysical sce- 
nario. 

In practice, we have considered the evolution of a nonrotat- 
ing neutron star with a gravitational mass of 2.75M , which 
is chosen to sit on the unstable branch of the equilibrium con- 
figurations and is endowed with an initial poloidal magnetic 
field of strength B c = 5 x 10 14 G extending also in the ex- 
terior space. As for the previous stellar solutions, we use a 
poly tropic EOS with T = 2 and K = 364.7 for the initial 
data and continue to use this isentropic EOS also for the sub- 
sequent evolution. The evolutions have been carried out in 
a computational domain with outer boundary at R out = 241 
km and a resolution of Ax = 0.11 km, corresponding to 272 
points covering the finest grid which extends up to ±15 km. 

Because the magnetic energy is only a small fraction of the 
binding energy, the hydrodynamical and spacetime evolution 
of the fluid star as it collapses to a BH is very similar to the 
unmagnetized case and this has been discussed in great detail 
in D9ll . The most important difference, therefore, is in the 
dynamics of the magnetic field, and this is shown in Fig. [T2l 
which reports two-dimensional cuts on the (x, z) plane of the 
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FIG. 12: Two-dimensional cuts on the (x, z) plane of the collapse to a BH of a magnetized NS. Shown with colors are the rest-mass density 
(colorcode from white to red) and the radial poynting vector (colorcode from blue to green) in units of 10 34 , while thin lines reproduce the 
magnetic-field lines. The different snapshots refer to times t = 0, 0.32, 0.57, 0.65, 1.0 and 1.1 ms, and an apparent horizon is marked with 
a thin red line starting from t — 0.57 ms. Note that all the matter is accreted into the hole and that a quadrupolar QNM ringdown is clearly 
visible in the Poynting flux. 
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FIG. 13: The same as the three bottom panels of Fig.[l2]but with a linear scale of 15 km to highlight the dynamics near the horizon. It is now 
very clear that a closed set of magnetic field lines is built just outside the horizon at t = 1.0 ms, that is radiated away as QNM of the BH. 



collapse to a BH of a magnetized NS. Shown with colors are 
the rest-mass density (colorcode from white to red) and the 
radial poynting vector (colorcode from blue to green), while 
thin solid lines reproduce the magnetic-field lines. 

At early times the star remains close to its initial state with 
the exception of a small transient induced by truncation error, 
which produces a small radiative outburst at t < 0.3 ms. As 



the instability to gravitational collapse develops, there is a re- 
arrangement of the external electromagnetic fields, driven by 
a toroidal electric field w —v v Bq produced in the inte- 
rior of the perfectly conducting star, and which is continuous 
across the stellar surface. As the collapse proceeds, the rest- 
mass density in the center and the curvature of the spacetime 
increase until an appararent horizon is found at t = 0.57 ms 
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FIG. 14: The same as in Fig.[l2j but where in addition to the rest-mass density (colorcode from white to red) and the magnetic-field lines (thin 
solid lines) we show the electrically-dominated regions (i.e. B 2 — E 2 < 0, colorcode from light blue to white in units of 10 23 ). 



and is marked with a thin red line in Fig. [12] (we have used the 
apparent-horizon finder described in [60]). 

As the stellar matter is accreted onto the BH (the rest-mass 
outside the horizon Mb, ut = is zero by t > 0.62 ms), 
the external magnetic field which was anchored on the stel- 
lar surface becomes disconnected, forming closed magnetic- 
field loops which carry away the electromagnetic energy in the 
form of dipolar radiation. This process, which has been de- 
scribed through a simplified non-relativistic analytical model 
in Ref. ll52ll . predicts the presence of regions where \E\ > \B\ 
as the toroidal electric field propagates outwards as a wave. 
This process can be observed very clearly in Fig. [131 which 
displays the same three bottom panels of Fig. [12] on a smaller 
scale of only 15 km to highlight the dynamics near the hori- 
zon. In particular, it is now very clear that a closed set of 
magnetic field lines is built just outside the horizon at t = 1.0 
ms, that is radiated away. Note also that our choice of gauges 
(which are the same used in Il6lll ) allows us to model without 
problems also the solution inside the apparent horizon. While 
the left panel of Fig. [13] shows that most of the rest-mass is dis- 
sipated away already by t = 0.65 ms (see discussion in ll62ll 
about why this happens), some of the matter remains on the 
grid near the singularity, anchoring there the magnetic field 
which slowly evolves as shown in the middle and right pan- 
els. A complementary view of the collapse process is also 
offered by Fig.QS which reports, in addition to the rest-mass 
density (colorcode from white to red) and the magnetic-field 




FIG. 15: Top Panel: Luminosity calculated at a distance r — 89 
km from the compact object. The black dotted line represents the 
time at which the apparent horizon is formed and the black dashed 
line corresponds to the time at which all the matter is well within 
the horizon. Bottom Panel: Evolution of the total radiated energy 
normalized to the initial magnetic energy. 
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FIG. 16: Left Panel: QNM ringdown of the magnetic field as measured through the magnetic flux at r = 37 km. Again, the black dotted line 
represents the time at which the apparent horizon is formed and the black dashed line corresponds to the time at which all the matter is well 
within the horizon; the dot-dashed line represents instead our fit to an exponential decay. Right Panel: Logarithm of the absolute values of the 
magnetic and electric fluxes as normalized to the initial magnetic flux. 



lines (thin solid lines), also the electrically-dominated regions 
(i.e. B 2 — E 2 < 0, colorcode from light blue to white). The 
larger scales used in this case makes it easier to follow the dy- 
namics of the closed field lines that once produced near the 
horizon, propagate as dipolar radiation at infinity. 

The total electromagnetic luminosity L ra d emitted during 
the collapse and computed as surface integral of the Poynt- 
ing flux over a spherical surface at 89 km is shown in the top 
panel of Fig. Q31 Note the presence of a rise during the col- 
lapse and of several pulses after the stellar matter has been 
accreted onto the black hole. The vertical dotted line repre- 
sents the time at which the apparent horizon is first found, 
while the vertical dashed line corresponds to the time at which 
all the matter is within the horizon (i.e. Mb, ut = 0). The 
peaks in the electromagnetic luminosity correspond to the 
closed magnetic-field loops that disconnect from the star and 
trasport electromagnetic energy. The bottom panel of Fig.[T5l 
on the other hand, reports the evolution of the total electro- 
magnetic energy lost in radiation E ra ^ and when normalized 
to the value of the initial magnetic energy outside the star, 
Eo • Our results indicate therefore a total electromagnetic effi- 
ciency which is ~ 5%; this result is smaller than the estimate 
made in Ref. ll5lll (which was of c± 20%), but, besides the dif- 
ferent initial data used, this difference can be easily accounted 
for by the fact that the gravitational collapse simulated here is 
considerably slower (and hence inefficient) than the OS one 
computed in Dill , where matter is free falling. Our efficiency 
is also smaller than the one computed in Ref. [52] and which 
is ~ 16% once the same definition for Eq is used. However, 
many other factors could be behind this difference, e.g., differ- 
ences in the initial data (use of a dipole everywhere in contrast 



to a dipole only outside the star as in our case), differences 
in the stellar models, differences in the numerical approach 
(treatment of the surface of the star of the transition between 
ideal and force-free MHD). A closer comparison between the 
two approaches will be carried out in a separate work. 

After BH formation, the luminosity decreases exponen- 
tially in a fashion which is typical of the QNM ringing of an 
electrovacuum electromagnetic field in a Schwarzschild BH 
spacetime. These QNMs are clearly visible also in the (ab- 
solute value of the) magnetic flux shown in the left panel of 
Fig. [16J from which a comparison with the perturb ative ex- 
pectations can be made. More specifically, by fitting the har- 
monic oscillations of the ringdown and the exponential de- 
cay we have computed the frequencies of the "ringing-down" 
magnetic-field flux for the I = 1 mode to be uj = 0.344054 — 
i 6.46731 kHz, corresponding to a nonrotating black hole of 
2.74 Mq. The agreement with the analytical value is excel- 
lent, with a relative error of only ~ 0.7% for the real part of 
the frequency and ~ 5.6% for the imaginary one 116311 . 

Finally, as a measure of the accuracy of our simulation we 
can compare the magnetic flux with the corresponding electric 
flux, which should vanish in the continuum limit since no net 
electric charge should be present. This is indeed the case, as 
can be deduced from the right panel of Fig. [T6J which reports 
the two fluxes normalized to the initial magnetic flux. Note 
that the electric flux is about 30 orders of magnitude smaller 
than the magnetic flux before BH formation, increasing after 
an apparent horizon is found, but remaining 15-10 orders of 
magnitude smaller. 
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V. CONCLUSIONS 

We have introduced a general-relativistic resistive MHD 
formalism as an extension of the special relativistic resistive 
MHD formalism reported in Ref. fl9ll for a 3+1 decompo- 
sition of the spacetime. Our numerical implementation has 
been made within the Cactus computational infrastructure 
as a continuation of the already existing general-relativistic 
hydrodynamics code Whisky |56|,[59|] and of the ideal-MHD 
code WhiskyMHD fTTIl . 

Our numerical approach exploits Implicit-Explicit (IMEX) 
methods and allows us to treat astrophysical problems in 
which different spatial regions fall into different regimes of 
conductivities. The flexibility introduced by using the Runge- 
Kutta will allow us to consider not only more general Ohm 
laws and a variety of astrophysical dynamos I24l64ll . but also 
to use better dispersion relations to calculate the velocities 
in the HLLE method and to describe more accurately non- 
relativistic systems ll65ll . 

Our implementation has been tested for a number of strin- 
gent tests and its robustness has been verified. The special- 
relativistic tests involved the propagation of circularly polar- 
ized Alfven waves, the evolution of current sheets and shock- 
tubes in one dimension, cylindrical and spherical explosion 
tests in two and three dimensions respectively, the evolution 
of stable and the collapse of unstable magnetized stars in dy- 
namical spacetime. We have compared our numerical results 
either with the analytical solution (in the cases where one ex- 
ists), or with the numerical ideal-MHD solution (in the limit 
of high conductivity), proving that our implementation is suit- 
able to describe regions with a wide range of conductivities, 
with or without large discontinuities and shocks. 

We have also considered genuinely general-relativistic tests 
in terms of the evolution of nonrotating magnetized stars ei- 
ther with fixed or fully dynamical spacetimes. Our stars have 
been endowed with magnetic fields of varying strength, ei- 
ther confined in their interior or permeating also the exterior 
space, and have been modelled with a non-uniform conductiv- 
ity that allows us to recover the ideal-MHD limit in the interior 
of the star and such the electrovacuum limit outside the star. 
All of our results indicate that the resistive implementation is 
able to follow the evolution of the oscillations triggered by the 
small truncation errors and that the associated eigenfrequen- 
cies match well those either reported with other hydrodynam- 
ics and ideal-MHD codes (HI EH] or from perturbation theory. 



Finally, we have considered the challenging and compre- 
hensive test represented by the gravitational collapse of a mag- 
netized nonrotating star to a BH. This scenario has an as- 
trophysical interest of its own as it could lead to the emis- 
sion of electromagnetic radiation, potentially detectable. In- 
deed we have found that as the collapse proceeds, electrically 
dominated regions develop and lead to the development of 
magnetic-field loops that propagate at the speed of light, car- 
rying away electromagnetic energy. Up to 5% of the initial 
magnetic energy can be lost in this way and the following 
evolution of the magnetic field follows a clean exponential 
decay, as expected by an electromagnetic perturbation in a 
Schwarzschild spacetime. The match of the measured QNMs 
and the perturbative predictions is well of a few percent or 
less. 

Our new code is now ready to be applied to study a vari- 
ety of astrophysical scenarios. These include the modeling of 
the magnetosphere that could be produced after the merger of 
binary neutron stars, or when the hypermassive neutron star 
collapses to a BH and is surrounded by a hot torus. The work 
in Ref. lf30h has already reported that under these conditions 
strong magnetic fields can be produced and that a jet-like mag- 
netic structure can develop. It is exciting to consider whether 
the resistive losses that are expected in the process will pro- 
vide sufficient energy to launch of a powerful jet, not yet ob- 
served in Ref. [30]. Also of great interest is to study BH mag- 
netospheres and the origin of jets so as to answer the ques- 
tion of whether an ergosphere is critical for the development 
of the Blandford-Znajek mechanism. Finally, our approach is 
also well suited to study the properties of accretion disk onto 
BHs and to elucidate the role that resistive losses play on the 
whole energetic budget. We will report on these applications 
in forthcoming works. 
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